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Introduction 


Abstract 


Calinaga (Moore 1857) is a rare and enigmatic Asian butterfly genus whose 
phylogenetic placement within Nymphalidae has only recently been established. The 
evolutionary history of Calinaga species however remains unknown. Here we explore 
the phylogeography of Calinaga using 1310 bp of sequence data from two molecular 
(mtDNA barcode and ribosomal protein S5 nuclear gene) and two morphological traits 
(genitalia and wing pattern). Within the proposed phylogenetic framework, we estimate 
the ages of divergence within the genus and reconstruct their historical biogeography. We 
found strong support for monophyly of Calinaga and support for the most recent accepted 
Species in the genus. Our results indicate that the common ancestor of Calinaga first 
split in the Eocene (~43 million years ago) in southern China, probably as a consequence 
of geological and environmental impacts of the collision of the Indian and Asian 
subcontinents. In the Oligocene/Miocene, the extrusion of Indochina from the continent 
caused further dramatic orogenetic changes that promoted isolation and speciation events 
within the genus while Pleistocene climatic changes also influenced the distribution and 
further speciation. A dispersal—vicariance analysis suggests that vicariance events have 
played a far more important role than dispersal in the distribution of extant species. 


under “Strip III., with Chilopodiform or Scolopendriform 
Larvae” alongside other nymphalids, with C. buddha as 


Calinaga (Moore, 1857) is an enigmatic butterfly genus 
mainly distributed in the Indochina region, but also ex- 
tending to Taiwan and the Himalayas (Fig. 1A). It is a 
prime example of a pattern that is found in diverse Hi- 
malayan-Indochinese genera of butterflies (e.g. Bhuta- 
nitis, Byasa, Callerebia, Loxerebia, Hemadara, Neope, 
Chrysozephyrus) and in many other groups including ver- 
tebrates (e.g. birds - Phylloscopus: Martens et al. 2011; 
frogs - tribe Paini: Che et al. 2010; plants - Roscoea: Zhao 
et al. 2015). Species of Calinaga inhabit forested areas at 
altitudes between about 1000-2000m where their larvae 
feed on various species of Morus (Lee 1958, Ashizawa 
and Muroya 1967, Chou 1994, Zhang and Hu 2012). 

The systematics of the genus has long been a topic for 
taxonomic discussion. Moore (1857) described the genus 


its type species. Kirby (1871) assigned Calinaga spe- 
cies to the Papilionidae (genus Parnassius), but he later 
(1877) placed the genus in the Nymphalidae (Hypolimnas 
group). Moore (1895) included Calinaga in the monoba- 
sic subfamily Calinaginae that is currently comprised of 
only this genus. Later Fruhstorfer (in Seitz 1927) recog- 
nized the genus under Nymphalidae based on certain mor- 
phological characters (genitalia, antennae, wing patterns), 
while Clark (1947) placed it in the “family” Danaidae and 
later gave it a separate family status, Calinagidae, between 
the “families” Argynnidae and Danaidae (Clark 1948). 
Ehrlich (1958a) recognized it as a distinct monotypic 
subfamily within the Nymphalidae, and, on the basis of 
characters of the adult integumental anatomy, suggested a 
close relationship to the subfamilies Satyrinae and Morph- 
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inae. This affinity was subsequently supported by studies 
on larval stages (Lee 1958) and adult morphology (Freitas 
and Brown 2004), as well as by recent molecular and phy- 
logenetic studies that have firmly placed the subfamily as 
the sister to Satyrinae + Charaxinae (Wahlberg and Wheat 
2008, Pefia and Wahlberg 2008, Wahlberg et al. 2009, 
Heikkila et al. 2012). At species level, Ehrlich (1958b) 
suggested that Calinaga was monotypic, comprised of a 
single species, C. buddha (Moore, 1857). The main prob- 
lem of the traditional taxonomic species-level framework 
for Calinaga, from Fruhstofer (in Seitz 1927) to Ehrlich 
(1958b), has been the practice of placing the vast major- 
ity of available names as subspecies or synonymous un- 
der C. buddha. Presently, taxonomists consider the genus 
composed from 7 to 11 species (D’Abrera 1985, 1993, 
Beccaloni et al. 2003, Lang 2012, Savela 2015, Sondhi et 
al. 2016, Vane-Wright, personal communication). A mo- 
lecular study to shed light on systematic and evolutionary 
history of these species has never been attempted. 

Located in the continental portion of Southeast Asia, ly- 
ing to the southeast of the Himalayas and south of China, 
the Indochina bioregion has long been recognized as an 
area with globally important levels of biodiversity (Myers 
et al. 2000). Recent meta-analyses of geological, climat- 
ic, and biological datasets have demonstrated impressive- 
ly that it is one of two major evolutionary hotspots for 
Southeast Asian biodiversity for a diverse range of plants 
and vertebrates (de Bruyn et al. 2014). It is a geologi- 
cally and topographically complex region, with intricate 
current and historical climatic patterns that contribute to 
its rich biotic diversity (e.g. Fontaine and Workman 1997, 
Hall 1998, 2001, An 2000, Morley 2000, Sterling et al. 
2006). Although several studies have examined the pro- 
cesses that have led to present patterns of species richness 
in plants and vertebrates (Bain and Hurley 2011, de Bruyn 
et al. 2014), little is known about the underlying causes for 
this remarkably high diversity, mechanisms of speciation, 
or origin of the Indochinese insect fauna. Until recently, 
some contributions addressing these topics were focused 
on species groups with high dispersal ability, such as but- 
terflies (Monastyrskii and Holloway 2013). 

For the first time in this study, we conduct an analysis of 
molecular data (mitochondrial COI-5’ and ribosomal pro- 
tein S5 nuclear gene (RpS5)) and we examine morphologi- 
cal characters (genitalia and wing pattern) for 51 specimens 
of Calinaga from widely distributed sites in the Indochina 
region (Fig. 1A, Table 1), for inferring preliminary infor- 
mation on the molecular phylogeny and the biogeographic 
processes which determined present-day distribution of the 
genus. An upcoming study by the first author (VT) will in- 
vestigate the type material held at NHM London, with the 
aim of adding new information to the existing dataset. 

Given the scarcity of information on the evolutionary 
history of Calinaga, the geographic and taxonomic cov- 
erage in our datasets will help to advance understanding 
of the diversification of this charismatic insect group in 
the Indochina bioregion and of other genera with similar 
distribution patterns. 
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Materials and Methods 


We examined 51 Calinaga individuals from 29 localities 
across India, South China (Yunnan, Shaanxi, Qinghai, Si- 
chuan), Laos, Vietnam, Myanmar and Thailand (Fig. 1A, 
Table 1). Photographs of all specimens used in the study 
as Well as detailed collecting data are available in the Bar- 
code of Life Data System (BOLD) (Ratnasingham and 
Hebert 2007) at https://doi.org/10.5883/DS-CALI. The 
specimens were initially identified as including represen- 
tatives of six species: C. aborica (N= 2), C. buddha (N= 
25), C. buphonas (N= 2), C. davidis (N= 6), C. lhatso (N= 
2) and C. sudassana (N= 14). Samples were submitted 
for sequencing using the original identifications, but they 
were re-identified following collection of genetic and 
morphological data (see paragraph 2.1). We were unable 
to obtain sequences from, among others, the following 
taxa: C. buddha buddha (Northern India), C. gautama 
(Sikkim), C. davidis davidis, C. fokienensis, C. kunagtun- 
gensis and C. pacifica (China), C. formosana (Taiwan), 
and C. distans (Southern Vietnam). 


Morphology 


Morphological identifications in this study were made 
based on illustrations of type material in Moore (1857), 
Melvill (1893), Leech (1892), de Niceville (1897), Seitz 
(1909, 1927), Tytler (1915), Oberthtr (1920), Lang 
(2012) and Sondhi et al. (2016). Male genitalia were dis- 
sected, prepared and photographed using traditional tech- 
niques (Supplementary Material S1: Sbordoni V. personal 
communication). 


Amplification and sequencing analysis 


DNA was extracted from one leg of each individual on 
Biomek FX liquid handling robot using a semi-automated 
DNA extraction protocol (Ivanova et al. 2006) on glass 
fiber plates (PALL Acroprep 96 with 3um GF membrane 
over 0.2um Bioinert membrane). DNA was eluted in 35- 
40ul of ddH,0 pre-warmed to 56°C and stored at -80°C. 
Polymerase chain reactions (PCR) and DNA sequencing 
were carried out following standard procedures for Lep- 
idoptera (Hajibabaei et al. 2005, deWaard et al. 2008). A 
658-bp fragment of cytochrome c oxidase subunit I (COI) 
was targeted for amplification using the primers LepF (3’ 
ATTCAACCAATCATAAAGATATTGG3’) and LepR 
(5° TAAACTTCTGGATGTCCAAAAAATCA3’). 

The entire nuclear RpS5 was initially amplified using 
the primers RpS5f/r (Wahlberg and Wheat 2008) and af- 
terwards a fragment covering 483 bp was obtained by nest- 
ed PCR using RpS5f and a novel specific primer (5’AC- 
MGTDCGYCGTCARGCTGTRGAYGTBTCWCC3’), 
developed from conserved regions. The PCR conditions 
were carried out as in Todisco et al. (2010). 

Sequences were obtained by using an ABI 3730x!I se- 
quencer following the manufacturer’s recommendations 
and they were edited and assembled using CodonCode 
Aligner 6.0.2. All sequences of Calinaga were submitted 
to GenBank (Accession Numbers in Table 1). 
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Figure 1. (A) Approximate geographic distributions (Shir6zu 1960, Lang 2012) and sampling localities (circles) for the species 
of Calinaga included in this study (with the exception of the sample CBUD-INDIN for which we do not have an exact locality). 
Species as initially identified are highlighted and shown in different colours. Note that many of these initially attributed names 
subsequently proved erroneous. The map was obtained using Quantum GIS 2.8.2 based on a map from Natural Earth (www.natu- 
ralearthdata.com). (B) Median-Joining Network of mtDNA. Circle size proportional to haplotype frequency; number of nucleotide 
substitutions indicated along connections, except for single or double substitutions. In both figures the species are highlighted and 
shown in different colours as initially identified. 


Phylogenetic analyses and fossil calibration for eight outgroups for species belonging to closely relat- 


Two sequences from C. buddha (COIL: EF683684, — ed lineages of Nymphalidae (Charaxinae and Satyrinae) 
AY090208; RpS5: EU141406) and one from C. davidis from Wahlberg and Wheat (2008) were retrieved from 
(COL: HQ658143), as well as COI and RpS5 sequences GenBank and added to our dataset (see Table 1). 
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Table 1. Taxon= Calinaga species identified in this study. Sample numbers (N), collection locality, population code, mtDNA hap- 


lotype code and GenBank Accession numbers. 


Population | Haplotype GenBank Accession 
taxon N Collection locality ede code COI KX233- RpS5 KX283- 
aborica 2 |Burma: N Sagaing, Tailing Hka river CABO-BURSA lal 590/612 390/398 
IHats6 1 China: Shaanxi, Tsnling Mts, Foping nature CLHA-CHSTE 570 379 

reserve, 1600 m 
davidis 2 |China: Qinghai, Ningyou Jiuzhi country CDAV-CHNJQ 576/581 382/384 
davidis nubilosa 2 |China: Qinghai, Ningyou Jiuzhi country CBUD-CHQNJ BITISOS 383/391 
davidis 1 |China: Sichuan, Yingxiuwan, 1000 m CDAV-CHSY| 596 394 
rr. China: W Sichuan, Dayi district, 23Km north 
davidis nubilosa 1 of Huashuiwan, 1500m CBUD-CHSDH nea 603 - 
AU nines l China: W Sichuan, Shimian—Mianning road, CBUD-CHSSM 4 599 _ 
2300 m 20 
Gavidisiiubidsa » | cie\e ine: Feehan, SOkmv- ef, Manian, CBUD-CHSTI H 580 hs 
1250 m ai 
davidis nubilosa 3 |China: W Sichuan, Yaan—-ergaz Shan, 1400 m| CBUD-CHSYS | H,,,H,, | 601/614/609 - 
Vitis 5 ul S Sichuan, SW of Dechang Miyi, 2100 CBUD-CHSDM we 7 
davidis buphonas | 3. [enina: Yunnan, Wer 1Okm SE, 2000-2800 | CDAV-CHYWE 613/584/616 387/399 
davidis buphonas 1 China: Yunnan, Surroundings of Dali city, CBUD-CHYDA 504 392 
2400 m 
os China: Yunnan, Ningjing Shan Tse—Kou 
davidis 2 Mekong River, 2100 m CBUD-CHYNM H 619/605 401/395 
davidis genestieri at ge PUMP AEN ENINE Wes SUAnaG SU see00 CBUD-CHYNS H 583 386 
davidis genestieri 1 |China: Yunnan, Wei-Xi, 2200 m CBUD-CHYWX elas 618 - 
davidis 1 a Yunnan, 40km N of Xhinxrong, 1850 CBUD-CHYXH Ho 586 0 
sudassana go) nina, MnpaneseWecn aleve SMA | eCAUECHYSA | eH 611/587 397 
Fugong, 1300 m ? 
davidis Sa aia ulna aN See e ar ON Dla || ESBUG_CHYSea|) Bae 617/620 400/402 
Gongshan, 1675 m 17 
oa China: N Yunnan, Sanba: Blotopel 
davidis cercyon 2 (Zhongdian), 2400 m CBUP-CHYSB It 598/591 = 
davidis buptionas | <2) |enina’NYunnan,, Zam NW of Tyang, 2000") -CBUD:CHYLI 5 578/571 E 
ihatso l China: N Yunnan, 55km N of Zhongdian, CLHA-CHYZH 592 7 
2450 m z 
sudassana 1 |NE India CBUD-INDIN 579 - 
sudassana 2 |India: Nagaland Naga hills CBUD-INDNNH 608/597 ~ 
sudassana 1 |Laos: Samneua CSUD-LAOSAM ye 607 - 
sudassana 1 | Thailand 602 - 
sudassana 1 | Thailand: Wiang Papao, Chiang Rai 574 - 
sudassana 2 |N Thailand: Samoeng CSUD-TAISAM 604/610 - 
sudassana SB ee a ote ol Meo ay aCeOISU ICL. gle CSUD-VNDVM 589/573/575 | 389/380/381 
Glang, 1700 m 
Vietnam: Hoang Lien national Park Sapa, Cat 595/606/588 
suyeoale : Cat village and Muong Hoa river, 1350 m a 600/585/582 eps eae 
davidis 1 |China EF683684 - 
davidis 1 |Stratford Butterfly Farm, UK AY090208 EU141406 
China: Tianma National Nature Reserve in 
eset 1 | Jinzhai County, Anhui ee 7 
eee 1 “ EU141357 FU141390 
wakefieldi 
Caligo telamonius | 1 - AY090209 EU141414 
Anaea troglodyta 1 - DQ338573 EU141428 
AICRACORIERONaS | 4 z AY090220 FU141424 
demophon 
Bicyclus anynana 1 - AY218238 EU141374 
Melanitis leda 1 - AY090207 EU141408 
Morpho peleides 1 - AY090210 EU141407 
Charaxes castor 1 - AY090219 EU141422 
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We used NETWORK 5 (Bandelt et al. 1999) to cal- 
culate a Median Joining (MJ) network representing the 
genealogical relationships among mtDNA haplotypes 
(Fig. 1B). Genetic distances were calculated in MEGA 
6 under the Kimura 2-parameter model of base substitu- 
tion (Tamura et al. 2013). The program JModelTest 2.1.7 
(Guindon and Gascuel 2003, Darriba et al. 2012) was 
employed to evaluate models of nucleotide substitutions. 
Only the best-fit models were subsequently used for esti- 
mations of tree topology. 

We used the BEAST package 1.8.3 (Drummond et al. 
2012) for Bayesian phylogeny reconstruction and diver- 
gence time estimation for concatenated genes. The anal- 
ysis was allowed to run for 50 million generations and 
the first 5000 trees were discarded as burnin. A consen- 
sus tree was calculated with posterior probability limit of 
0.5, maximum clade credibility tree and median heights. 
Monophyly was enforced on a) Satyrinae, b) Charaxinae, 
and c) Satyrinae + Charaxinae (Fig. 2). The tree prior 
was set to Yule Process Speciation model with a random 
starting tree. The analysis was performed using an un- 
correlated relaxed clock (lognormal relaxed distribution 
not selected) for both genes. Trees were visualized using 
FIGTREE 1.4 (Rambaut 2009). 

Recent phylogenetic studies calibrated with fossil data 
(Pefia and Wahlberg 2008; Wahlberg et al. 2009) suggest 
that Nymphalidae began to diversify in the late Creta- 
ceous around 90 Mya, satyrines (Calinaginae + Satyrinae 
+ Charaxinae) around 75 Mya and Satyrinae + Charax- 
inae about 70 Mya. Thus, in our BEAST analyses the tree 
root was calibrated with normal distribution at 75 Mya 
(stdev 3), Satyrinae + Charaxinae at 70 Mya (stdev 3.5) 
and Charaxes + Euxanthe at 22 Mya (stdev 1) (see Fig. 2). 


Tests of demographic equilibrium and population ex- 
pansion in mtDNA haplogroup 


Haplotype and nucleotide diversity were calculated us- 
ing DnaSP 5.0 (Librado and Rozas 2009). Demograph- 
ic equilibrium was tested in each mtDNA haplogroup 
by calculating F, (Fu 1997) and R, (Ramos-Onsins and 
Rozas 2002) statistics, which have been shown to provide 
sensitive signals of population expansion (Ramos-Onsins 
and Rozas 2002). Arlequin 3.0 (Excoffier et al. 2005) and 
DnaSP 5.0 were used to compute F, and R,, respectively, 
and to test their statistical significance by simulating ran- 
dom samples (10,000 replicates) under the null hypoth- 
esis of selective neutrality and constant population size 
using coalescent algorithms (both modified from Hudson 
1990). P-values for the two statistics were obtained as the 
proportion of simulated values smaller than or equal to 
the observed values (a= 0.05). Expected mismatch dis- 
tributions and parameters of sudden expansion t= 2 u 
T were calculated using Arlequin 3.0 by a generalized 
least-squares approach (Schneider and Excoffier 1999), 
under models of pure demographic expansion and spatial 
expansion (Ray et al. 2003, Excoffier 2004). The prob- 
ability of the data under the given model was assessed 
by the goodness-of-fit test implemented 1n Arlequin 3.0. 
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Parameter confidence limits were calculated in Arlequin 
3.0 through a parametric bootstrap (1000 simulated ran- 
dom samples). 


Ancestral Area Reconstruction 


Dispersal-vicariance optimization of ancestral areas im- 
plemented in DIVA 1.1 (Ronquist 1996) was used to infer 
the biogeographic history of the group. The distributions 
of populations were divided into five biogeographic units 
designated following Che et al. (2010) and Proches and 
Ramdhani (2012), as follows: A) Southwestern China 
ecozone, B) Himalaya-Tibetan plateau region, C) North- 
ern Sino-Himalaya, D) Southern Sino-Himalaya, E) In- 
dochina (Fig. 2). The least number of dispersal events 
(12) was achieved when maxareas was set to 3. 


Results 


Morphology 


The morphological study on the wing pattern of our spec- 
imens, based on illustrations in the most relevant litera- 
ture, suggested the presence of four species: C. aborica, 
C. davidis, C. lhatso and C. sudassana (Table 1). Many of 
the samples procured as C. buddha proved to be misidenti- 
fications of either C. davidis or C. sudassana. In addition, 
the samples of C. buddha (COI. EF683684, AY090208; 
RpS5: EU141406) and C. davidis (COI: HQ658143) re- 
trieved from GenBank were re-identified respectively as 
C. davidis and (possibly) C. /actoris (Table 1). 

Male genitalia in Calinaga were found to be very uni- 
form and thus genitalic characters were of little use in 
discrimination of species. Although Lang (2012) empha- 
sized the shape and ‘pointiness’ of the dorsal apex of the 
valve as an important character, we did not find this trait 
useful or consistent within species. All the adult and gen- 
italia images are available in Supplementary Material S1. 


Phylogenetic analyses and network analysis 


For COI (658 bp), 54 sequences (51 sequences analyzed 
in this study plus 3 sequences retrieved from GenBank) 
showed 23 different haplotypes (Fig. 1B; Table 1) with 71 
variable sites. The 24 RpS5 sequences (483 bp) displayed 
4 haplotypes and 6 variable sites. 

Bayesian phylogeny analysis was used to reconstruct 
phylogenetic relationships (Fig. 2) for the concatenat- 
ed genes. The resulting tree was based on an alignment 
of 1310 bp and it was rooted using outgroup sequences 
of members of two subfamilies: Charaxinae (Charaxes 
castor, Anaea troglodyte, Archaeoprepona demophon, 
Euxanthe wakefieldi) and Satyrinae (Caligo telamonius, 
Morpho peleides, Bicyclus anynana, Melanitis leda) (Ta- 
ble 1). The selected model of evolution for RpSS nuclear 
gene was GTR+G+I (Rodriguez et al. 1990) with estimat- 
ed base frequencies, 4 gamma categories, partitioned into 
positions (1+2),3 and for COI was HKY (Hasegawa et al. 
1985) with estimated base frequencies, no heterogeneity 
model and not partitioned. 
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Figure 2. Bayesian phylogeny of Calinaga estimated in BEAST using concatenated data. Purple squares are calibration points (root: 75 
+ 3; Satyrinae + Charaxinae 70 + 3.5, Charaxes + Euxanthe 22 + 1). Monophyly was enforced on nodes marked with orange squares. 
The inset map shows the biogeographic regions used in DIVA analysis: A) Southwestern China ecozone, B) Himalaya-Tibetan plateau 
region, C) Northern Sino-Himalaya, D) Southern Sino-Himalaya, E) Indochina. Colored dots correspond to haplogroups on the tree. 


The Bayesian tree for concatenated genes resulted in 
a well-supported Calinaga clade (Fig. 2). Beside distinct 
lineages for C. [hatso (N= 1) and C. aborica (N= 2), two 
major clades were observed. The first clade, consisting of 
specimens identified as C. buddha (N= 20), C. davidis (N= 
6) and C. buphonas (N= 2), was called Haplogroup A “da- 
vidis”’ because all samples were re-identified as C. davidis 
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(Table 1). The latter also included two GenBank sequenc- 
es for C. buddha (COI: EF683684, AY090208; RpSS: 
EU141406). Examination of the voucher images for these 
sequences showed that these also were indeed C. davidis 
(Min et al. 2008; Pefia and Malm 2012). The second clade, 
consisted of specimens initially identified as C. sudassana 
(N= 14), C. buddha (N= 5) and C. [hatso (N= 1) and it was 
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called Haplogroup B “sudassana’” because all samples 
were subsequently re-identified as C. sudassana (Table 1). 

The MJ network analysis (Fig. 1.B) of the mtDNA 
haplotypes was fully consistent with our Bayesian phylo- 
genetic analysis. Forty-five mutational steps joined Hap- 
logroup A “davidis” (including all C. davidis sequences; 
Fig. 1B) and Haplogroup B “sudassana” (including all 
C. sudassana sequences with C. /hatso sequences CL- 
HA-CHSFT1). The genetic differentiation between the 
two haplogroups, estimated as the mean distance (Kimura 
2 parameter uncorrected p-distance) between groups, was 
10.3%. The analysis also highlighted a star-like config- 
uration of Haplogroup A “davidis” with the H,. haplo- 
type shared among Yunnan populations localized in the 
headwater area for the three major rivers flowing from 
Tibet (Salween, Mekong and Yangtse; Fig. 1A, Fig. 2). 
In this geographical area Callinga displays its highest 
genetic variability, and Haplogroup A “davidis” and B 
“sudassana’ coexist (Fig. 2). Thirty-four and thirty-three 
mutational steps joined C. aborica with Haplogroup A 
and B respectively. In addition, in Haplogroup A “da- 
vidis”, a GenBank sequence for C. davidis (HQ658143) 
from Tianma National Nature Reserve in Jinzhai County, 
Anhui, China (Xia and Hao 2011), possibly representing 
the taxon C. Jactoris Fruhstorfer 1908, was divergent 
(Fig. 1B). Finally, in the Haplogroup B “sudassana’”, 
the single specimen CLHA-CHSFT1 from Tsinling Mts, 
Foping nature reserve in Shaanxi, initially identified as C. 
lhatso based on morphological similarity, showed identi- 
cal COI and RpS5 sequences to C. sudassana. Extraction 
and amplification of this sample was conducted twice in 
two different laboratories and generated the same results. 
Eleven mutational steps (Fig. 1B) separated the DNA bar- 
code sequence of this individual from the other C. /hat- 
so from Yunnan (CLHA-CHYZH1) despite their similar 
wing patterns. 


Tests of demographic equilibrium and estimation of 
divergence times 


According to Fs and R, statistics, calculated for Hap- 
logroup A “davidis” (N= 20) and Haplogroup B “sudas- 
sana” (N= 30) mtDNA sequence sets (Table 2), the null 
hypothesis of constant population size could be rejected 
for these two phylogeographic units (bold in Table 2). 
Mismatch distribution of these groups was examined ac- 
cording to sudden expansion model (Table 2), and good- 
ness of fit tests did not show significant deviations from 
expected distributions, so that parameter T= 2ut could be 
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used to estimate the time (t) elapsed from population ex- 
pansion: estimated values of t and their 5% and 95% con- 
fidence limits are shown in Table 2. According to Brower 
(1994) proposed mutation rates in Lepidoptera (u= 0.01 
substitutions/site/lineage/Ma), demographic expansion of 
the Haplogroup B “sudassana” could be traced back to t= 
260 kya (126-403 kya). A similar estimate was obtained 
for Haplogroup A “davidis” t= 168 kya (70-362 kya). 

The most recent common ancestor of Calinaga has 
been estimated to have split from Satyrinae + Charaxinae 
sometime between 67-84 Mya (Wahlberg et al. 2009). 
Results of our BEAST analysis reveal that the first major 
split within Calinaga occurred in the Eocene, around 43 
Mya (63-25 Mya) (Fig. 2). The time of the most recent 
ancestor ((MRCA) of C. aborica + “davidis group” was 
estimated to have occurred in the Oligocene ~30 Mya 
(47-17 Ma), whereas it was ~21 Mya (38-8 Ma) for C. 
lhatso + “sudassana group”. Finally, the tMRCA for both 
davidis and sudassana groups was estimated at ~12 Mya 
(24-4 Ma). 


Biogeographic ancestral reconstructions 


Results of the ancestral area reconstructions showed that 
the common ancestor of Calinaga probably occurred 
before ~43 Mya in area ADE (Fig. 2), corresponding to 
Southwestern China ecozone, Southern Sino-Himalaya 
and Indochina. Before ~30 Mya this also represented the 
ancestral area for clade C. aborica + Haplogroup A “da- 
vidis”. Subsequently in the Miocene before ~15 Mya, in 
a more restricted area ranging between the Southwestern 
China ecozone and Southern Sino-Himalaya (AD), the 
common ancestor of Haplogroup A “davidis” was pres- 
ent with subsequent dispersal (in the late Miocene and 
Pleistocene) into the Himalaya-Tibetan Plateau region. 
The common ancestor of clade C. /hatso + Haplogroup B 
“sudassana’ colonized the Southern Sino-Himalaya and 
Indochina (AE) before ~20 Mya. Finally, before ~10 Mya 
Indochina (E) represented the ancestral area of the Hap- 
logroup B “sudassana’. 


Discussion 


Taxonomic implications 


Despite six species initially identified (Fig. 1, Table 1), 
our molecular results identified only four main groups 
(C. aborica, C. sudassana, C. davidis and C. lhatso). 
These were also consistent with subsequent careful mor- 


Table 2. N, number of mtDNA sequences; H, number of unique mtDNA haplotypes; h, mtDNA haplotype diversity (4SD); a, 
mtDNA nucleotide diversity (4SD). Tests of demographic equilibrium and mismatch analysis in mtDNA phylogeographic groups: 
FS, Fu’s FS statistic; R,, Rozas and Ramon-Onsins’ R, statistic; t, sudden demographic expansion parameter, with 5% and 95% 
confidence limits; t, true time since population expansion, from t = 2 u T (where p is the substitution rate per gene). Significantly 


small values of F, and ——— are in bold. 


Haplogroup | N | 


_ X10) | Fs | R, | | r%) [7 (95%)| tka) | t6%) | t (95%) 


Loss|, £6: an 260 126 403 


30 tt toe 885+0.042 | 0.45940.575 | -4.203 | 0.121 | 3.867 
B 20} 7 | 0.840.068 0.288+40.3 | -1.369 | 0.143 | 2.537 


1.061 | 5.498 168 70 S63 
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phological identifications (Fig. 2, Table 1). This is quite 
reflective of the persistent confusion about the taxonomic 
status and number of species within Calinaga. Histori- 
cally, various authors have recognized different numbers 
and sets of valid Calinaga species, from one (Ehrlich 
1958b) or three (Fruhstorfer 1914) to seven (D’Abrera 
1985, 1993, Beccaloni et al. 2003), ten (Savela 2015) or 
eleven species (Sondhi et al. 2016). Many of the samples 
we procured as C. buddha proved to be misidentifications 
of either C. davidis or C. sudassana. A thorough examina- 
tion of the type material of all the names under Calinaga 
is required to clarify whether some of these names should 
be synonymized. Within our examined specimens we 
found that C. /hatso and C. aborica were morphological- 
ly and genetically distinct entities. However, the identity 
of the specimen identified as C. /hatso (CLHA-CHSFT1), 
sharing its haplotype with C. sudassana, remains unclear. 
Many of the specimens initially identified to sub-species 
under C. davidis showed little or no genetic differentia- 
tion (e.g. C. davidis nubilosa, C. davidis genestieri, C. 
davidis cercyon, C. davidis buphonas etc.). 


Historical biogeography 


The question of where and when Calinaga originated and 
how its species are related have not been previously in- 
vestigated. Although no fossils are known for the genus, 
previous phylogenetic studies on Nymphalidae that have 
used fossils for molecular calibration have dated the last 
common ancestor of Calinaga + (Satyrinae + Charax- 
inae) at 67-84 Mya (Pefia and Wahlberg 2008; Wahlberg 
et al. 2009). Since these studies used only a single Calin- 
aga species, the age of the crown group of Calinaga and 
divergences within the genus have remained unknown. 
Our results show that the first split in Calinaga occurred 
in Eocene at around 43 Mya (Fig. 2). The large gap be- 
tween this date and the age of the most recent common 
ancestor of Calinaga and Satyrinae + Charaxinae (67-84 
Mya) suggests diversification after the Cretaceous/Tertia- 
ry boundary with subsequent extinction of sister lineages 
(Ree and Smith 2008). The initial collision between India 
and Asia in the early Cenozoic (50-55 Ma) set in motion 
major orogenic associated environmental changes that 
continued through the Oligocene, Miocene and into mod- 
ern times which were likely a major driving force in evo- 
lution of Calinaga, and numerous other plant (e.g. Liu et 
al. 2002, Mao et al. 2010) and animal lineages (e.g. Ruber 
et al. 2004, Zhang et al. 2006, Nazari et al. 2007, Che et 
al. 2010, Zhao et al. 2015) found in the Himalayan-Ti- 
betan Plateau. In fact, the evolutionary history of organ- 
isms in this region helps elucidate parallel geological and 
climatic developments, and often aids the recognition of 
unrecorded geological events. 

Our results suggest that the common ancestor of Cali- 
naga probably occurred before ~43 Mya (63-25 Mya) in 
the region now corresponding to the Southwestern China 
ecozone, Southern Sino-Himalaya and Indochina (ADE: 
Fig. 2). Since the Indochinese peninsula did not exist be- 
fore ~22 Mya (Leloup et al. 2001), the ancestral distribu- 
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tion of Calinaga was likely confined to the Southwestern 
China ecozone (D) and Southern Sino-Himalaya (A). 
Since its center of genetic diversity seems to lie on the 
south-eastern border of the Tibetan Plateau, in particular 
where the three major rivers sourced in Tibet (Salween, 
Mekong and Yangtse rivers; Fig. 1A, Fig. 2) are in prox- 
imity, Calinaga may well have originated in the Tibetan 
region. Geological records indicate that around the Oligo- 
cene/Miocene boundary (~23 Mya), the rapid extrusion 
of Indochina and the cumulative Cenozoic deformation 
of South-East Asia caused by underthrusting of the In- 
dian Plate beneath the Eurasian Plate, resulted in intense 
orogeny in South-East Tibet (Zhang et al. 2013). In par- 
ticular in the Early Miocene (~20 Mya), driven by rapid 
deformation of the eastern Himalayan syntaxis (Robinson 
et al. 2014), large modifications were induced in the river 
systems of this region (e.g. Salween, Mekong, Yangtze 
and Yarlung Tsangpo-Irrawaddy river: Zhang and Hu 
2012). These geological upheavals have likely provided 
ample opportunities for allopatric isolation and subse- 
quent speciation in Calinaga. 

Our molecular dating suggests that C. /hatso and Hap- 
logroup B “sudassana” split into distinct lineages in Early 
Miocene ~21 Mya (38-8 Mya), possibly isolating the last 
common ancestor of C. /hatso in the Southern Sino-H1- 
malaya (A). Southern Sino-Himalaya and Indochina (AE) 
possibly represented the ancestral area for C. /hatso + Hap- 
logroup B “sudassana’. Considering that before ~22 Mya 
the Indochinese peninsula did not exist, it is more likely 
that the ancestral area of this clade was mainly in Southern 
Sino-Himalaya (A). The common ancestor of Haplogroup 
B “sudassana” began diversifying around ~12 Mya (24-4 
Mya) and it probably expanded its range from Southern 
Sino-Himalayan region (A) to Indochina peninsula (E). 
Conversely, the origin of C. aborica + Haplogroup A “da- 
vidis” must be traced back to Oligocene (~30 Mya), when 
the first uplift of the Himalayan—Tibetan Plateau occurred 
(29 — 65 Mya: Zhao et al., 2015). Afterwards, in Miocene 
~12 Mya (21-4 Mya), Haplogroup A “davidis” diversified 
and expanded its range into Central-East China. 

In addition, our analyses show that Haplogroup A “da- 
vidis” and Haplogroup B “sudassana” experienced demo- 
graphic expansion in the Pleistocene, ~168 kya (70-363 
kya) and ~260 ka (16-403 kya) respectively. Therefore, the 
genetic diversity and the present-day distributions of Cali- 
naga should be attributed not only to the dramatic geolog- 
ical events in Oligocene/Miocene, but also to the impacts 
of the cycles of glacial advance/retreat that occurred in this 
region during the Pleistocene (Quan et al. 2014). 

Finally, the taxon C. formosana, for which we exam- 
ined several specimens but could not obtain any sequenc- 
es (Supplementary Material S1), is morphologically al- 
lied to the Haplogroup A “davidis”. Since Taiwan formed 
~4-5 Mya at a complex convergent boundary between the 
Philippine Sea Plate and the Eurasian Plate (“The Penglai 
Orogeny”: Teng 1990, Liu et al. 2000), the most plausible 
hypothesis is that C. formosana must have dispersed from 
mainland China to Taiwan recently in Pleistocene. 
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Conclusion 


The genus Calinaga probably originated in the South- 
East Tibet in Eocene following the immense geological 
and environmental impact caused by the collision be- 
tween Indian and Asian subcontinents. The extrusion 
of Indochina from the continent during the Oligocene/ 
Miocene further prompted dramatic orogenetic changes 
that promoted isolation and speciation events in the ge- 
nus. More recently, in the Pleistocene, climatic changes 
further modified the distribution of species and probably 
facilitated vicariant speciation events. 

Since we did not sample or sequence specimens from 
all of the available names under Calinaga, we cannot 
make any definitive statements about the number of valid 
Species warranted to be recognized as such, although the 
existence of many superfluous names is evident. From the 
names of the genus and the species coined by early British 
lepidopterists including F. Moore, it is apparent that they 
drew inspiration from Hindu mythological characters. In 
Sanskrit, Naga refers to mythical reptilian creatures found 
in Indian religions (Hinduism, Buddhism and Janism) who 
were often worshipped as deities. Among them, “Kaliya” 
(or Kalya, “Kalia-Naga”, Calinaga) was a particularly no- 
torious and poisonous one living in Yamuna river in Vrin- 
davan (Uttar Pradesh). After an encounter with Krishna, 
Kaliya surrendered and was sent to exile (Bhagavata Pu- 
rana, 16:10). It seems that the modern taxonomy of Cali- 
naga is in need of a Krishna to conquer these superfluous 
names and cleanse its taxonomy albeit after careful exam- 
ination of the types and sequencing of additional material. 
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